In this vignette we will introduce how to fit Whittle–Matérn fields
with general smoothness based on finite element and rational
approximations. The theory for this approach is provided in Bolin et al. (2022) and Xiong,
Simas, and Bolin (2022). For the implementation, make use of
the rSPDE package for the rational approximations.
These models are thus implemented using finite element approximations. Such approximations are not needed in for integer smoothness parameters, and for the details about the exact models we refer to the vignettes
For further details on the construction of metric graphs, see Working with metric graphs
We begin by loading the rSPDE and
MetricGraph packages:
As an example, we consider the following metric graph
library(sp)
line1 <- Line(rbind(c(0,0),c(1,0)))
line2 <- Line(rbind(c(0,0),c(0,1)))
line3 <- Line(rbind(c(0,1),c(-1,1)))
theta <- seq(from=pi,to=3*pi/2,length.out = 20)
line4 <- Line(cbind(sin(theta),1+ cos(theta)))
Lines = SpatialLines(list(Lines(list(line1),ID="1"),
Lines(list(line2),ID="2"),
Lines(list(line3),ID="3"),
Lines(list(line4),ID="4")))
graph <- metric_graph$new(lines = Lines)
graph$plot()To construct a FEM approximation of a Whittle–Matérn field with general smoothness, we must first construct a mesh on the graph.
In the command build_mesh, the argument h
decides the largest spacing between nodes in the mesh. As can be seen in
the plot, the mesh is very coarse, so let’s reduce the value of
h and rebuild the mesh:
The next step is to build the mass and stiffness matrices for the FEM basis.
We are now ready to specify the model \[
(\kappa^2 - \Delta)^{\alpha} \tau u = \mathcal{W}
\] for the Whittle–Matérn field \(u\). For this, we use the
matern.operators function from the rSPDE
package:
sigma <- 1.3
r <- 0.2 # range
nu <- 0.9
rspde.order <- 2
kappa <- sqrt(8 * nu) / r
op <- matern.operators(C= graph$mesh$C, G = graph$mesh$G, d = 1,
nu = nu, kappa = kappa, sigma = sigma,
m = rspde.order)As can be seen in the code, we specify \(\kappa\) via the practical correlation
range \(\sqrt{8\nu}/\kappa\). Also, the
model is not parametrized by \(\tau,
\alpha\) but instead by \(\sigma,
\nu\). Here, sigma denotes the standard deviation of
the field and nu is the smoothness parameter, which is
related to \(\alpha\) via the relation
\(\alpha = \nu + 1/2\). The object
op_cov contains the matrices needed for evaluating the
distribution of the stochastic weights in the FEM approximation.
Let us simulate the field \(u\) at the mesh locations and plot the result:
If we want to evaluate \(u(s)\) at
some locations \(s_1,\ldots, s_n\), we
need to multiply the weights with the FEM basis functions \(\varphi_i(s)\) evaluated at the locations.
For this, we can construct the observation matrix \(\boldsymbol{\mathrm{A}}\), with elements
\(A_{ij} = \varphi_j(s_i)\), which
links the FEM basis functions to the locations. This can be done by the
function mesh_A in the metric graph object. To illustrate
this, let us simulate some observation locations on the graph and
construct the matrix:
obs.per.edge <- 50
obs.loc <- NULL
for(i in 1:graph$nE) {
obs.loc <- rbind(obs.loc,
cbind(rep(i,obs.per.edge), runif(obs.per.edge)))
}
n.obs <- obs.per.edge*graph$nE
A <- graph$mesh_A(obs.loc)In the code, we generate \(50\) observation locations per edge in the graph, drawn at random. It can be noted that we assume that the observation locations are given in the format \((e, d)\) where \(e\) denotes the edge of the observation and \(d\) is the position on the edge, i.e., the relative distance from the first vertex of the edge.
To compute the precision matrix from the covariance-based rational
approximation one can use the precision() method on object
returned by the matern.operators() function:
As an illustration of the model, let us compute the covariance function between the process at the mid point of the second edge and all other locations in the mesh. The covariances can be calculated as \[ \overline{\boldsymbol{\mathrm{A}}} \boldsymbol{\mathrm{Q}}^{-1}\overline{\boldsymbol{\mathrm{v}}}. \] Here, \(\boldsymbol{\mathrm{Q}}\) is the precision matrix obtained from the rational approximation, \(\boldsymbol{\mathrm{A}}\) is an identity matrix since we are evaluating the approximation in the nodes of the FEM mesh, \(\overline{\boldsymbol{\mathrm{v}}}\) is the \((m+1)\)-fold vertical concatenation of the vector \(\boldsymbol{\mathrm{v}}\), where \(\boldsymbol{\mathrm{v}}\) is a vector with all basis functions evaluated in \(s=0.5\).
There is built-in support for computing log-likelihood functions and
performing kriging prediction in the rSPDE package which we
can use for the graph model. To illustrate this, we use the simulation
to create some noisy observations of the process. We generate the
observations as \(Y_i = u(s_i) +
\varepsilon_i\), where \(\varepsilon_i
\sim N(0,\sigma_e^2)\) is Gaussian measurement noise.
Let us now fit the model. To this end we first must compute the loglikelihood function as function of the parameters we want to estimate. We define the loglikelihood function parametrized using the logarithm of each parameter to avoid constrained optimization.
We now set some suitable initial values for the optimization and fit
the model using optim.
theta0 <- graph_starting_values(graph, range_par = TRUE,
nu = TRUE, like_format = TRUE,
manual_data = Y, log = TRUE)
theta <- optim(theta0, mlik, method = "L-BFGS-B")
print(data.frame(sigma = c(sigma, exp(theta$par[1])),
range = c(r, exp(theta$par[2])),
nu = c(nu, exp(theta$par[3])),
sigma.e = c(sigma.e, exp(theta$par[4])),
row.names = c("Truth", "Estimates")))## sigma range nu sigma.e
## Truth 1.300000 0.200000 0.9000000 0.1000000
## Estimates 1.385231 0.243925 0.8771482 0.1047377
Given that we have estimated the parameters, let us compute the kriging predictor of the field given the observations at the mesh nodes.
Let us update the model object with the fitted parameters:
sigma_est <- exp(theta$par[1])
kappa_est <- exp(theta$par[2])
nu_est <- exp(theta$par[3])
op_cov <- update(op,
user_kappa = kappa_est,
user_sigma = sigma_est,
user_nu = nu_est)We can now perform kriging with the predict()
method:
Since we predicted at the mesh nodes, Aprd was chosen as
the identity matrix with the size of the mesh. If we wish to predict at
some other locations, we can simply change Aprd to be the
observation matrix for those locations.
The estimate is shown in the following figure
Let us now illustrate how to simulate a data set with replicates and
then fit a model to such data. To simulate a latent model with
replicates, all we do is set the nsim argument to the
number of replicates.
Now, let us generate the observed values \(Y\):
Note that \(Y\) is a matrix with 20
columns, each column containing one replicate. Now, the remaining of the
code is identical to the previous case. The
rSPDE.matern.loglike() function automatically identifies
the replicates from the fact that \(Y\)
is a matrix with more than one column.
theta0 <- graph_starting_values(graph, range_par = TRUE,
nu = TRUE, like_format = TRUE,
manual_data = Y.rep, log = TRUE)
mlik <- rSPDE.construct.matern.loglike(op_cov, Y=Y.rep, A=A, parameterization = "matern")
theta <- optim(theta0, mlik, method = "L-BFGS-B")
print(data.frame(sigma = c(sigma, exp(theta$par[1])),
range = c(r, exp(theta$par[2])),
nu = c(nu, exp(theta$par[3])),
sigma.e = c(sigma.e, exp(theta$par[4])),
row.names = c("Truth", "Estimates")))## sigma range nu sigma.e
## Truth 1.300000 0.2000000 0.9000000 0.3000000
## Estimates 1.275934 0.1898274 0.8898979 0.3052531
We also have an R-INLA implementation of the rational
SPDE approach for metric graphs.
We begin by defining the model by using the
rspde.metric_graph() function. This function contains the
same arguments as the function rspde.matern(). We refer the
reader to the R-INLA
implementation of the rational SPDE approach vignette for further
details.
We begin by adding the observations to the graph:
graph$add_observations(data = data.frame(y=Y,
edge_number = obs.loc[,1],
distance_on_edge = obs.loc[,2]),
normalized = TRUE)Let us create the model object:
By default, the order of the rational approximation is 2.
We can now create the A matrix and the index object by
using the same functions as for the rational SPDE approach, namely,
rspde.make.A() and rspde.make.index(),
supplying the graph object as the mesh argument. Here we create the
A matrix:
Now, let us create the index object:
The remaining is standard: we create the formula object, the stack
object, and then fit the model by using the inla()
function. So, first we create the formula object:
Now we create the inla.stack object:
stk.dat <- inla.stack(
data = graph_data_spde(rspde_model), A = Abar, tag = "est",
effects =
c(
rspde.index,
list(Intercept = 1)
)
)Finally, we can fit the model:
rspde_fit <- inla(f.s, data = inla.stack.data(stk.dat),
control.inla = list(int.strategy = "eb"),
control.predictor = list(A = inla.stack.A(stk.dat), compute = TRUE),
inla.mode = "experimental"
)We can use the same functions as the rspde fitted models
in inla. For instance, we can see the results in the
original scale by creating the result object:
## mean sd 0.025quant 0.5quant 0.975quant mode
## std.dev 1.408070 0.1800640 1.091420 1.395020 1.797600 1.367280
## range 0.241322 0.0524461 0.153785 0.236118 0.358954 0.225903
## nu 0.902068 0.1309470 0.654682 0.898850 1.165900 0.889826
Let us compare with the true values:
result_df <- data.frame(
parameter = c("std.dev", "range", "nu"),
true = c(sigma, r, nu),
mean = c(
result_fit$summary.std.dev$mean,
result_fit$summary.range$mean,
result_fit$summary.nu$mean
),
mode = c(
result_fit$summary.std.dev$mode,
result_fit$summary.range$mode,
result_fit$summary.nu$mode
)
)
print(result_df)## parameter true mean mode
## 1 std.dev 1.3 1.4080696 1.3672847
## 2 range 0.2 0.2413216 0.2259025
## 3 nu 0.9 0.9020683 0.8898255
We can also plot the posterior marginal densities with the help of
the gg_df() function:
posterior_df_fit <- gg_df(result_fit)
library(ggplot2)
ggplot(posterior_df_fit) + geom_line(aes(x = x, y = y)) +
facet_wrap(~parameter, scales = "free") + labs(y = "Density")R-INLA implementationWe will do kriging on the mesh locations:
Let us now add the observations for prediction:
graph$add_observations(data = data.frame(y=rep(NA,nrow(pred_loc)),
edge_number = pred_loc[,1],
distance_on_edge = pred_loc[,2]),
normalized = TRUE)Let us compute the A matrix at the prediction
locations:
Let us build the prediction stack and gather it with the estimation
stack. To this end, we set the argument only_pred to
TRUE, in which it will return the data.frame
containing the NA data.
ef.prd <-
c(rspde.index, list(Intercept = 1))
stk.prd <- inla.stack(
data = graph_data_spde(rspde_model, only_pred=TRUE),
A = Abar_prd, tag = "prd",
effects = ef.prd
)
stk.all <- inla.stack(stk.dat, stk.prd)Let us obtain the predictions:
rspde_fitprd <- inla(f.s,
data = inla.stack.data(stk.all),
control.predictor = list(
A = inla.stack.A(stk.all),
compute = TRUE, link = 1
),
control.compute = list(
return.marginals = FALSE,
return.marginals.predictor = FALSE
),
control.inla = list(int.strategy = "eb")
)Let us now extract the indices of the predicted nodes and store the means:
id.prd <- inla.stack.index(stk.all, "prd")$data
m.prd <- rspde_fitprd$summary.fitted.values$mean[id.prd]Finally, let us plot the predicted values. To this end we will use
the plot_function() graph method.
R-INLA implementation to fit models with
replicatesLet us begin by cloning the graph and clearing the observations on the cloned graph:
We will now add the data with replicates to the graph:
graph_rep$add_observations(data = data.frame(y=as.vector(Y.rep),
edge_number = rep(obs.loc[,1], n.rep),
distance_on_edge = rep(obs.loc[,2], n.rep),
repl = rep(1:n.rep, each = n.obs)),
group = "repl",
normalized = TRUE)To fit the model with replicates we need to create the corresponding
A matrix and index object:
Abar.rep <- rspde.make.A(
mesh = graph_rep, index = rep(1:n.obs, times = n.rep),
repl = rep(1:n.rep, each = n.obs)
)
mesh.index.rep <- rspde.make.index(
name = "field", mesh = graph_rep,
n.repl = n.rep
)Similarly, let us create a new rspde model object:
Let us now create the corresponding inla.stack object,
where we set the repl argument in the function
graph_data_spde to __all since we want to use
all replicates:
st.dat.rep <- inla.stack(
data = graph_data_spde(rspde_model_rep,
repl = "__all"),
A = Abar.rep,
effects = mesh.index.rep
)Observe that we need the response variable y to be a
vector. We can now create the formula object, remembering
that since we gave the name argument field, when creating
the index, we need to pass field.repl to the
formula:
We can, finally, fit the model:
rspde_fit_rep <-
inla(f.rep,
data = inla.stack.data(st.dat.rep),
family = "gaussian",
control.predictor =
list(A = inla.stack.A(st.dat.rep)),
inla.mode = "experimental"
)We can obtain the estimates in the original scale with the
rspde.result() function:
## mean sd 0.025quant 0.5quant 0.975quant mode
## std.dev 1.272090 0.0270196 1.220920 1.271370 1.32702 1.269320
## range 0.187300 0.0104929 0.169241 0.186351 0.21024 0.183601
## nu 0.904973 0.0686166 0.760356 0.909706 1.02751 0.926811
Let us compare with the true values of the parameters:
result_rep_df <- data.frame(
parameter = c("std.dev", "range", "nu"),
true = c(sigma, r, nu),
mean = c(
result_fit_rep$summary.std.dev$mean,
result_fit_rep$summary.range$mean,
result_fit_rep$summary.nu$mean
),
mode = c(
result_fit_rep$summary.std.dev$mode,
result_fit_rep$summary.range$mode,
result_fit_rep$summary.nu$mode
)
)
print(result_rep_df)## parameter true mean mode
## 1 std.dev 1.3 1.2720948 1.2693174
## 2 range 0.2 0.1873000 0.1836008
## 3 nu 0.9 0.9049733 0.9268108
We can also plot the posterior marginal densities with the help of
the gg_df() function:
posterior_df_fit_rep <- gg_df(result_fit_rep)
ggplot(posterior_df_fit_rep) + geom_line(aes(x = x, y = y)) +
facet_wrap(~parameter, scales = "free") + labs(y = "Density")inlabru implementationThe inlabru package allows us to fit models and do
kriging in a straighforward manner, without having to handle
A matrices, indices nor inla.stack objects.
Therefore, we suggest the reader to use this implementation when using
our implementation to fit real data.
Let us clear the graph, since it contains NA
observations we used for prediction, add the observations again, and
create a new rSPDE model object:
graph$clear_observations()
graph$add_observations(data = data.frame(y=Y,
edge_number = obs.loc[,1],
distance_on_edge = obs.loc[,2]),
normalized = TRUE)
rspde_model <- rspde.metric_graph(graph)Let us now load the inlabru package and create the
component (which is inlabru’s formula-like object). Since
we are using the data from the graph, inlabru will also
obtain the locations from the graph, thus, there is no need to provide
the locations. However, we need a name for the locations for using
inlabru’s predict method. Therefore, we can choose any name
for the location that is not a name being used in the graph’s data. In
our case we will use the name loc:
Now, we can directly fit the model, where we pass the name of the
location variable as the loc argument of the
graph_data_spde() function:
rspde_bru_fit <-
bru(cmp,
data=graph_data_spde(rspde_model, loc = "loc"),
options=list(
family = "gaussian",
inla.mode = "experimental")
)Let us now obtain the estimates of the parameters in the original
scale by using the rspde.result() function:
## mean sd 0.025quant 0.5quant 0.975quant mode
## std.dev 1.408250 0.1799740 1.091840 1.395180 1.797670 1.367240
## range 0.241471 0.0526009 0.153772 0.236222 0.359531 0.225915
## nu 0.901841 0.1315030 0.653310 0.898657 1.166680 0.889757
Let us compare with the true values of the parameters:
result_bru_df <- data.frame(
parameter = c("std.dev", "range", "nu"),
true = c(sigma, r, nu),
mean = c(
result_bru_fit$summary.std.dev$mean,
result_bru_fit$summary.range$mean,
result_bru_fit$summary.nu$mean
),
mode = c(
result_bru_fit$summary.std.dev$mode,
result_bru_fit$summary.range$mode,
result_bru_fit$summary.nu$mode
)
)
print(result_bru_df)## parameter true mean mode
## 1 std.dev 1.3 1.4082541 1.3672408
## 2 range 0.2 0.2414709 0.2259154
## 3 nu 0.9 0.9018408 0.8897565
We can also plot the posterior marginal densities with the help of
the gg_df() function:
posterior_df_bru_fit <- gg_df(result_bru_fit)
ggplot(posterior_df_bru_fit) + geom_line(aes(x = x, y = y)) +
facet_wrap(~parameter, scales = "free") + labs(y = "Density")inlabru implementationIt is very easy to do kriging with the inlabru
implementation. We simply need to provide the prediction locations to
the predict() method.
In this example we will use the mesh locations. To this end we will
use the get_mesh_locations() method. We also set
bru=TRUE and loc="loc" to obtain a data list
suitable to be used with inlabru.
Now, we can simply provide these locations to the
predict method along with the fitted object
rspde_bru_fit:
Finally, let us plot the predicted values. To this end we will use
the plot_function() graph method:
We can also use our inlabru implementation to fit models
with replicates. We will consider the same data that was generated
above, where the number of replicates is 30.
For this implementation we will use the rspde_model_rep
object.
We can now create the component, passing the vector with the indices
of the replicates as the replicate argument. For the
replicate argument we use the function
graph_repl_spde() to extract the vector of replicates,
where we set choose __all as the replicates, since we want
all replicates.
cmp_rep <-
y ~ -1 + Intercept(1) + field(loc, model = rspde_model_rep,
replicate = graph_repl_spde(rspde_model_rep, "__all"))Now, we are ready to fit the model. Observe that we set the
repl argument to __all since we will use all
replicates and the loc argument to loc, which
is the name we gave to the locations.
rspde_bru_fit_rep <-
bru(cmp_rep,
data=graph_data_spde(rspde_model_rep,
repl = "__all",
loc = "loc"),
options=list(
family = "gaussian",
inla.mode = "experimental")
)We can obtain the estimates in the original scale with the
rspde.result() function:
result_bru_fit_rep <- rspde.result(rspde_bru_fit_rep, "field", rspde_model)
summary(result_bru_fit_rep)## mean sd 0.025quant 0.5quant 0.975quant mode
## std.dev 1.280840 0.0273472 1.227870 1.280590 1.335290 1.280150
## range 0.191240 0.0106389 0.169550 0.191665 0.211052 0.193754
## nu 0.897091 0.0672550 0.787134 0.889103 1.046570 0.860258
Let us compare with the true values of the parameters:
result_bru_rep_df <- data.frame(
parameter = c("std.dev", "range", "nu"),
true = c(sigma, r, nu),
mean = c(
result_bru_fit_rep$summary.std.dev$mean,
result_bru_fit_rep$summary.range$mean,
result_bru_fit_rep$summary.nu$mean
),
mode = c(
result_bru_fit_rep$summary.std.dev$mode,
result_bru_fit_rep$summary.range$mode,
result_bru_fit_rep$summary.nu$mode
)
)
print(result_bru_rep_df)## parameter true mean mode
## 1 std.dev 1.3 1.2808406 1.2801461
## 2 range 0.2 0.1912403 0.1937536
## 3 nu 0.9 0.8970910 0.8602582
We can also plot the posterior marginal densities with the help of
the gg_df() function:
posterior_df_bru_fit_rep <- gg_df(result_bru_fit_rep)
ggplot(posterior_df_bru_fit_rep) + geom_line(aes(x = x, y = y)) +
facet_wrap(~parameter, scales = "free") + labs(y = "Density")